Small RNA and Degradome Sequencing in Floral Bud Reveal Roles of miRNAs in Dormancy Release of Chimonanthus praecox

Chimonanthus praecox (wintersweet) is highly valued ornamentally and economically. Floral bud dormancy is an important biological characteristic in the life cycle of wintersweet, and a certain period of chilling accumulation is necessary for breaking floral bud dormancy. Understanding the mechanism of floral bud dormancy release is essential for developing measures against the effects of global warming. miRNAs play important roles in low-temperature regulation of flower bud dormancy through mechanisms that are unclear. In this study, small RNA and degradome sequencing were performed for wintersweet floral buds in dormancy and break stages for the first time. Small RNA sequencing identified 862 known and 402 novel miRNAs; 23 differentially expressed miRNAs (10 known and 13 novel) were screened via comparative analysis of breaking and other dormant floral bud samples. Degradome sequencing identified 1707 target genes of 21 differentially expressed miRNAs. The annotations of the predicted target genes showed that these miRNAs were mainly involved in the regulation of phytohormone metabolism and signal transduction, epigenetic modification, transcription factors, amino acid metabolism, and stress response, etc., during the dormancy release of wintersweet floral buds. These data provide an important foundation for further research on the mechanism of floral bud dormancy in wintersweet.


Introduction
MicroRNAs (miRNAs), the key post-transcriptional regulators of eukaryotic gene expression, are a class of endogenous non-coding small RNAs (sRNAs) with 18-24 nucleotides (nt) in length [1]. Mature miRNAs are processed from larger primary miRNAs (pri-miRNAs) with a series of endonucleases and then associated with Argonaute (AGO) proteins to form RNA-induced silencing complexes (RISCs). Based on perfect or nearperfect complementarity between miRNAs and target mRNAs, RISCs negatively regulate gene expression by targeting mRNAs for cleavage or translational repression [2,3]. To date, a large number of miRNAs have been identified in various plant species, such as rice, maize, Arabidopsis, poplar, grape, and peach [4][5][6][7][8][9]. Furthermore, numerous studies have shown that miRNAs play important regulatory roles in various biological processes, including plant growth and development, signal transduction, and biotic and abiotic stress [10][11][12][13]. miRNAs show both conservation and diversity among plant lineages. Some miRNAs are conserved in angiosperms and even embryophyta, whereas more miRNAs show species specificity, which reflects the natures of rapid evolution and functional differentiation of miRNAs [14,15].
Bud dormancy is a temporary cessation of visible growth during plant growth and development, and it is a biological characteristic of perennial deciduous plants adapted 2 of 17 to environmental and seasonal changes during long-term evolution [16][17][18]. According to the factors that cause dormancy, bud dormancy is generally classified into three types: paradormancy, endodormancy, and ecodormancy [16]. Endodormancy, also known as natural dormancy, is controlled by endogenous factors of dormancy structure. Short days (SDs) and/or low temperatures are the main factors for inducing bud endodormancy, and endodormant buds do not resume growth until the chilling requirement is fulfilled [19]. In this study, endodormancy is referred to as dormancy unless specified otherwise.
Only several studies have identified miRNAs that regulate bud dormancy in trees, including pear, apple, peony, peach, and grape [21,[37][38][39][40]. In white pear, miRNA6390 directly targets PpDAM to inhibit its expression, thereby promoting the release of floral bud dormancy [21]. Furthermore, miR156k, miR159a, miR167a, etc., were found to be differentially expressed in tree peony buds between endodormancy and ecodormancy [38]. In peach, miR6285 was found to regulate dormancy release in flower buds by targeting asparagine-rich protein (NRP), which is involved in ABA signal transduction [39]. The miR159-MYB module could be involved in regulating apple bud dormancy by mediating the homeostasis of endogenous ABA [37]. Most of these data were obtained with high-throughput sequencing; thus, the regulation mechanisms of miRNAs in bud dormancy are unclear.
Chimonanthus praecox (wintersweet), a perennial woody ornamental plant, is a rare precious flower and cut flower material in winter that has high ornamental and economic value. In autumn, wintersweet floral buds are induced to enter dormancy, and a certain period of chilling accumulation is essential to break dormancy [41,42]. In recent years, the deficiency of chilling accumulation caused by global warming has interrupted normal dormancy release in some cases, which seriously affects the flowering process and quality of wintersweet. Although some studies on the molecular mechanisms of wintersweet floral bud dormancy regulation were performed at the transcriptomic and proteomic levels [41], the roles of miRNAs involved in chilling-mediated dormancy control of wintersweet floral buds remain unclear. High-throughput sequencing technology and bioinformatics are effective in identifying miRNAs. Therefore, in this study, sRNA sequencing (sRNA-seq) of wintersweet floral buds was performed to identify miRNAs involved in the regulation of dormancy release through next-generation sequencing. Meanwhile, the target genes of these miRNAs were identified with degradome sequencing (parallel analysis of RNA end sequencing, PARE-seq). This study provides a basis for further analyzing the roles of miRNAs during the dormancy of wintersweet floral buds and also provides a platform for investigation of the biological functions of miRNAs in wintersweet. In addition, our work contributes to a better understanding of the mechanism of bud dormancy in perennial plants.

Overview of sRNA Sequencing
With C. praecox 'Concolor' plants as the materials, six samples of floral buds (FBs) (FB.Nov, nF16, FB150, FB300, FB450, and IB570, three biological replicates for each sample) with different chilling accumulation during the dormancy process were used for sRNA sequencing, and 18 sRNA libraries were constructed. These libraries generated a total of 358,543,196 clean reads through high-throughput sequencing and preliminary filtering. After removing contamination, low-quality reads, and reads with length <18 nt, 324,509,314 clean tags were obtained, including 39,956,685 unique tags (non-redundant sequences). The sequence length of these sRNAs was concentrated at 20-24 nt, accounting for 86.92% of clean tags. For all samples, 21-and 24-nt sRNAs were the most abundant classes (Figure 1). These lengths are consistent with the specificity of Dicer-like (DCL1-4) proteins that cleave double-stranded RNA (dsRNA) to form sRNA [43], which is consistent with the length distribution of plant sRNAs. Among all samples, 24-nt sRNAs were the most abundant in FB.Nov, FB150, and FB450, whereas 21-nt sRNAs were the most abundant in the other three samples ( Figure 1). The results showed that the sRNAs' sequences were diverse and complex in these samples.

Overview of sRNA Sequencing
With C. praecox 'Concolor' plants as the materials, six samples of floral buds (FBs) (FB.Nov, nF16, FB150, FB300, FB450, and IB570, three biological replicates for each sample) with different chilling accumulation during the dormancy process were used for sRNA sequencing, and 18 sRNA libraries were constructed. These libraries generated a total of 358,543,196 clean reads through high-throughput sequencing and preliminary filtering. After removing contamination, low-quality reads, and reads with length <18 nt, 324,509,314 clean tags were obtained, including 39,956,685 unique tags (non-redundant sequences). The sequence length of these sRNAs was concentrated at 20-24 nt, accounting for 86.92% of clean tags. For all samples, 21-and 24-nt sRNAs were the most abundant classes (Figure 1). These lengths are consistent with the specificity of Dicer-like (DCL1-4) proteins that cleave double-stranded RNA (dsRNA) to form sRNA [43], which is consistent with the length distribution of plant sRNAs. Among all samples, 24-nt sRNAs were the most abundant in FB.Nov, FB150, and FB450, whereas 21-nt sRNAs were the most abundant in the other three samples (Figure 1). The results showed that the sRNAs' sequences were diverse and complex in these samples. These clean tags were mapped to the GenBank, Rfam, RepeatMasker, and C. praecox transcriptome databases successively for annotation and classification. A total of 11.45% and 22.87% of clean tags were annotated as non-coding structural RNAs (rRNA, tRNA, snRNA, and snoRNA) and transcriptome sequences (degraded fragments of wintersweet mRNAs), respectively ( Table 1). The annotated tags were removed, and the remaining clean tags were further analyzed to identify known and novel miRNAs of wintersweet.  These clean tags were mapped to the GenBank, Rfam, RepeatMasker, and C. praecox transcriptome databases successively for annotation and classification. A total of 11.45% and 22.87% of clean tags were annotated as non-coding structural RNAs (rRNA, tRNA, snRNA, and snoRNA) and transcriptome sequences (degraded fragments of wintersweet mRNAs), respectively ( Table 1). The annotated tags were removed, and the remaining clean tags were further analyzed to identify known and novel miRNAs of wintersweet.

Identification of Known and Novel miRNAs
All unannotated clean tags were aligned with known plant miRNAs sequences in miRBase 21.0. A total of 862 known miRNAs were identified in all six samples (18 libraries) of wintersweet (Table S1), containing 74,164,064 tags and accounting for 22.85% of total clean tags (Table 1). These known miRNAs belonged to 754 families (Table S2). The expression profiles of these known miRNAs showed great differences in the amounts of tags. Among them, the total abundance of miR166-y was the highest, with 31,649,892 tags from all libraries, followed by miR396-x (27,464,586). However, 79.70% of known miRNAs had fewer than 10 tags (Table S1). The results indicated that the expression of these known miRNAs could have spatiotemporal specificity, and that these known miRNAs might play different regulatory roles in the dormancy of wintersweet floral buds. In addition, only 27 precursors and their secondary structures of 19 miRNAs from 862 known miRNAs were identified in the C. praecox transcriptome database, which included the same precursors forming mature miRNAs with different sequences and mature miRNAs with the same sequences derived from different precursors (Table S3). The main reasons for this result were the lack of both an available genome database and a complete transcriptome database for wintersweet.
Further, the tags annotated as known miRNAs were excluded, and the remaining clean tags were aligned against the C. praecox transcriptome database to identify the loci that may serve as precursors (pre-miRNA) of miRNAs. By predicting and analyzing the hairpin structures of these putative precursors, 402 potential novel miRNAs were identified in all samples of wintersweet (Table S4), containing 1,014,923 tags and accounting for 0.31% of the total clean tags ( Table 1). Compared with known miRNAs, the expression levels of novel miRNAs were relatively low (Table S4). The precursors of all 402 novel miRNAs possessed typical stem-loop structures. These precursors were 57-371 nt in length; the minimum folding free energy (MFE) of these precursors hairpin structures ranged from −166.5 kcal·mol −1 to −18.0 kcal·mol −1 (Table S5), which was similar to the characteristic of other plant miRNAs precursors [44,45]. Among the precursors of 402 novel miRNAs, 39 precursors formed mature miRNAs on both 3 and 5 arms (Table S6). In addition, similar to known miRNAs, some mature novel miRNAs with the same sequence were derived from different precursors, such as novel-m0010 and novel-m0011 (Table S5). Secondary structure prediction of several precursors is shown in Figure S1.
Sequence analysis of all of the 1264 miRNAs (862 known and 402 novel) revealed that their lengths ranged from 18 nt to 26 nt. miRNAs with 21 nt were the most abundant, accounting for 85.50% of the total miRNAs, followed by 22-(7.12%) and 20-nt (7.00%) miRNAs, whereas only one 26-nt miRNA was identified ( Figure 2a). Studies have shown that the base composition of miRNAs can affect their physicochemical and biological properties [46]. When AGO proteins bind to miRNAs to form RISCs, different AGO proteins are preferentially enriched for the different initiating nucleotides of miRNAs [47].
The analysis of first nucleotide bias showed that the first nucleotide of 37.82% (478) and 30.06% (380) of miRNAs were U and A, respectively. For miRNAs of different lengths, the first nucleotides of 20-, 21-, and 22-nt miRNAs were biased toward U, whereas those of other miRNAs were biased toward A (Figure 2b).
OR PEER REVIEW 5 of 18 [47]. The analysis of first nucleotide bias showed that the first nucleotide of 37.82% (478) and 30.06% (380) of miRNAs were U and A, respectively. For miRNAs of different lengths, the first nucleotides of 20-, 21-, and 22-nt miRNAs were biased toward U, whereas those of other miRNAs were biased toward A ( Figure 2b).

Differential Expression of miRNAs during Dormancy Release
The expression levels of all identified miRNAs were calculated by using the transcripts per million (TPM) algorithm. According to the expression levels of the miRNAs in each sample, all samples were compared in pairs to screen the differentially expressed miRNAs (DEMs) between different samples. A total of 301 DEMs were obtained ( Figure  S2, Table S7). The results from the statistical analysis of DEMs between different samples are shown in Figure 3a, which reveals that the differences in miRNAs expression became more significant with chilling accumulation; notably, the differences between the dormancy breaking sample (IB570) and samples remaining in dormancy (FB.Nov, nF16, and FB150/300/450) were the greatest.

Differential Expression of miRNAs during Dormancy Release
The expression levels of all identified miRNAs were calculated by using the transcripts per million (TPM) algorithm. According to the expression levels of the miRNAs in each sample, all samples were compared in pairs to screen the differentially expressed miRNAs (DEMs) between different samples. A total of 301 DEMs were obtained ( Figure S2, Table S7). The results from the statistical analysis of DEMs between different samples are shown in Figure 3a, which reveals that the differences in miRNAs expression became more significant with chilling accumulation; notably, the differences between the dormancy breaking sample (IB570) and samples remaining in dormancy (FB.Nov, nF16, and FB150/300/450) were the greatest.
To identify miRNAs involved in the regulation of wintersweet floral bud dormancy release, we focused on the comparative analysis of DEMs among IB570/FB.Nov, IB570/nF16, IB570/FB150, IB570/FB300, and IB570/FB450, and 23 DEMs shared in the five comparisons were screened, including 10 known and 13 novel miRNAs (Figure 3b, Table S8). Of the 10 known miRNAs, miR157-x, miR167-x, miR167-y, miR390-y, and miR395-y were highly conserved and were studied in various plants; miR473 was reported only in several tree species (Populus trichocarpa, Prunus persica, and Citrus trifolia) [9,49,50]; the other miRNAs (miR2092, miR5179, miR6293, and miR9566) were rarely reported. A heatmap depicting the expression of the 23 DEMs showed that 15 miRNAs (5 known and 10 novel) were upregulated, and 8 miRNAs (5 known and 3 novel) were downregulated in IB570 compared with other samples (FB.Nov, nF16, and FB150/300/450) (Figure 3c). Among these DEMs, only miR167-x had very high expression (TPM > 20,000) in all samples, whereas the expression levels (TPM) of the remaining miRNAs were <100 (Table S8), indicating that miR167-x plays an important regulatory role in the induction and release of dormancy in wintersweet floral buds. miR6293-x and miR9566-y were barely expressed in all samples with remaining dormancy, except IB570 (Table S8), which indicated that miR6293-x and miR9566-y might promote dormancy release in wintersweet floral buds. With the continuous accumulation of chilling, the expression of novel-m0186-3p showed an increasing trend, whereas the expression of novel-m0077-5p showed a decreasing trend compared with that of FB.Nov (Table S8), indicating that novel-m0186-3p and novel-m0077-5p could respond to a chilling environment and have opposite roles in regulating the dormancy of wintersweet floral buds.
OR PEER REVIEW 6 of 18 novel-m0186-3p and novel-m0077-5p could respond to a chilling environment and have opposite roles in regulating the dormancy of wintersweet floral buds.

Summary of Degradome Sequencing
Equal amounts of total RNA from each sample were mixed for degradome sequencing, and one degradome library was constructed. After sequencing and data filtering were completed, 37,913,833 clean tags were obtained, including 12,290,544 unique tags. After removing non-coding structural RNAs and polyN fragments, the remaining clean tags were aligned against the C. praecox transcriptome database. A total of 31.33% of

Summary of Degradome Sequencing
Equal amounts of total RNA from each sample were mixed for degradome sequencing, and one degradome library was constructed. After sequencing and data filtering were completed, 37,913,833 clean tags were obtained, including 12,290,544 unique tags. After removing non-coding structural RNAs and polyN fragments, the remaining clean tags were aligned against the C. praecox transcriptome database. A total of 31.33% of the unique tags were positively aligned with the database (annotated as cDNA_sense), and 29.14% of the unique tags were reverse complemented to the database (annotated as cDNA_antisense) ( Table 2). Tags annotated as cDNA (cDNA_sense and cDNA_antisense) were used for the detection of miRNAs target genes.

Identification of miRNA Targets
Taking the C. praecox transcriptome database as reference, the cleavage sites of 1264 miRNAs identified with the described sRNA-seq were detected. A total of 220,567 cleavage sites of 1186 miRNAs were detected, corresponding to 33,310 degradation genes and 230,989 miRNAtarget pairs (Table S9). According to the abundance of tags in degradome sequencing, these miRNA-target pairs were divided into five categories: category 0, 1, 2, 3, and 4, which included 1374, 1312, 98,246, 8170, and 121,887 miRNA-target pairs, respectively (Figure 4a, Table S9). Based on the detection results of cleavage sites, target plots (T-plots) were drawn in units of genes. The T-plot of two target genes is shown in Figure 4c,d.
A total of 1532 of the 1707 target genes identified were annotated with corresponding functions. Among them, most encode metabolism-related enzymes or proteins that are involved in amino acid metabolism, epigenetic regulation (DNA methylation, histone modification, chromatin remodeling, and RNA modification, etc.), phytohormone metabolism and signal transduction (ABA, IAA, and ET, etc.), glucose metabolism, and stress response. Some target genes are related to flower development (such as FUL-like, SEP1, COL5, and FCA) and SPL, MYB, WRKY, and bHLH transcription factors. In addition, there are also a large number of target genes encoding structural component domains of certain proteins (Table S10). These results indicate that miRNAs are involved in the regulation of various biological processes during the dormancy release of wintersweet floral buds. Among the 1707 target genes, Gene Ontology (GO) enrichment analysis was performed on targets from miRNA-target pairs that were classified into categories 0, 1, and 2 with high reliability; unfortunately, only 320 of 855 genes were annotated. In Cellular component, significantly enriched GO terms were only found in the endomembrane system, microtubule associated complex, phosphatase complex, and protein serine/threonine phosphatase complex. The Molecular functions were mainly related to transition metalion binding, transferase activity,transferring glycosyl groups, isomerase activity, and anion binding. Regarding Biological process, these target genes were mainly involved in single-organism process, single-organism cellular process, single-organism metabolic process, biological regulation, small molecule metabolic process, cellular component organization or biogenesis, and cellular component organization ( Figure 5, Table S12). The GO enrichment analysis provided abundant resources for target gene isolation, localization, and expression profiling in wintersweet. Among the 1707 target genes, Gene Ontology (GO) enrichment analysis was performed on targets from miRNA-target pairs that were classified into categories 0, 1, and 2 with high reliability; unfortunately, only 320 of 855 genes were annotated. In Cellular component, significantly enriched GO terms were only found in the endomembrane system, microtubule associated complex, phosphatase complex, and protein serine/threonine phosphatase complex. The Molecular functions were mainly related to transition metalion binding, transferase activity, transferring glycosyl groups, isomerase activity, and anion binding. Regarding Biological process, these target genes were mainly involved in singleorganism process, single-organism cellular process, single-organism metabolic process, biological regulation, small molecule metabolic process, cellular component organization or biogenesis, and cellular component organization ( Figure 5, Table S12). The GO enrichment analysis provided abundant resources for target gene isolation, localization, and expression profiling in wintersweet.

Discussion
Wintersweet is a valuable flowering plant in winter. A certain period of chilling accumulation is necessary to break its floral buds' dormancy and promote flowering. miRNAs are important post-transcriptional regulators involved in various biological processes in plants. Here, sRNA sequencing and degradome sequencing were combined to study the roles of miRNAs in the chilling-mediated dormancy release of wintersweet floral buds for the first time.
In the present study, the transcriptome data of Chimonanthus praecox was used for the annotation of miRNAs and their target genes rather than the genome, though two genomes were published [52,53]. The main reasons for this include the following two aspects: the annotation information of both genomes that were deposited in the NCBI (National Center for Biotechnology Information) is lacking; moreover, there may be some conflicts of interest due to the copyright issues involved.

Discussion
Wintersweet is a valuable flowering plant in winter. A certain period of chilling accumulation is necessary to break its floral buds' dormancy and promote flowering. miRNAs are important post-transcriptional regulators involved in various biological processes in plants. Here, sRNA sequencing and degradome sequencing were combined to study the roles of miRNAs in the chilling-mediated dormancy release of wintersweet floral buds for the first time.
In the present study, the transcriptome data of Chimonanthus praecox was used for the annotation of miRNAs and their target genes rather than the genome, though two genomes were published [52,53]. The main reasons for this include the following two aspects: the annotation information of both genomes that were deposited in the NCBI (National Center for Biotechnology Information) is lacking; moreover, there may be some conflicts of interest due to the copyright issues involved.
In plants, 21-nt miRNAs are typical and play important roles in post-transcriptional gene silencing by mediating mRNA cleavage or translational repression [43]. Of the 1264 miRNAs (862 known and 402 novel) identified in this study, 21-nt miRNAs accounted for 85.50%, which was consistent with the report that most plant miRNAs are 21 nt long [60]. The analysis of first nucleotide bias showed that the first bases of 20-, 21-, and 22-nt miRNAs were biased toward U, and those of other miRNAs with different lengths were biased toward A (Figure 2). AGO proteins are the core components of RISCs and have certain bias toward the starting nucleotide of miRNAs when they associate with miRNAs to form RISCs. In Arabidopsis, AGO1 predominantly binds miRNAs that initiate with base U; AGO2 and AGO4 preferentially recruit miRNAs with 5 terminal A; whereas AGO5 is biased toward miRNAs with 5 terminal C [46,47,54]. These findings indicate that the nucleotide at the 5 terminal of a miRNA is closely related to its biological activity, and if it is changed, the miRNA may interact with one different AGO protein and show different biological functions.
In plants, nonconserved miRNAs are generally expressed at low levels or mainly expressed in specific cells or under specific growth conditions [7]. In this study, the abundance of novel miRNAs was generally lower than that of known miRNAs, which is consistent with the results of various plants reported so far, such as those from analyses of the Japanese pear (Pyrus pyrifolia "Kosui"), turnip rape (Brassica rapa L.), and hot pepper (Capsicum annuum L.) [45,46,59]. Furthermore, according to the miRBase database, these novel miRNAs have not been reported in other plants. Therefore, it is speculated that these novel miRNAs may be species-specific to wintersweet and are nonconserved miRNAs, whose expression requires certain induction conditions or have certain tissue specificity.
In the study, a comparative analysis of DEMs among IB570/FB.Nov, IB570/nF16, IB570/FB150, IB570/FB300, and IB570/FB450 was performed to identify miRNAs associated with the dormancy release of wintersweet floral buds, and 23 common DEMs were screened (Figure 3b). Among them, miR157-x, miR167-x, miR167-y, miR390-y, and miR395y were highly conserved and downregulated after dormancy release. During the whole dormancy process, miR167-x showed high expression, whereas the expression of miR167-y was relatively low, and the other three miRNAs were also expressed at low levels. These miRNAs were also identified in other plant studies of floral bud dormancy, such as for pear and peach plants; however, their expression characteristics are not always consistent among plants. In white pear, two miR157 (Pyr-miR157a and miR157) and Pyr-miR390f were upregulated after floral bud dormancy release, whereas the expression trends of two miR167 (Pyr-miR167g and Pyr-miR167j) and Pyr-miR395i showed no significant difference throughout the dormancy process; moreover, during dormancy, these miRNAs, excluding Pyr-miR395i, showed high expression levels [21]. In the floral buds of the Japanese pear, miR395 was not identified, and the expression levels of miR157, miR167, and miR390 did not show significant differences between endodormancy and ecodormancy; only miR167 showed high expression in both stages [59]. In addition, miR167b and miR167c of the peach were highly expressed in dormant floral buds, whereas their expression was at a lower level after dormancy release [39]. These results suggest that miRNAs from the same family have different expression characteristics in different plants. Some scholars have proposed that different plants have different responses to the same cold stress [61]. Which indicates that the roles of miRNAs and the regulation mechanisms of floral bud dormancy differ among plants in their long-term evolution.
Identifying the target genes and their functions is an important way to analyze the biological functions and mechanisms of miRNAs. In plants, miRNAs mainly silence gene expression at the post-transcriptional level by targeting mRNAs for cleavage [43]. In this study, we identified 1707 target genes of 21 DEMs in wintersweet floral buds via degradome sequencing, which included 1940 miRNA-target pairs.
Of the 1707 target genes, a fraction were involved in flower development. And among these genes, a FUL-like homologous gene, CpFUL-like, was targeted by novel-m0128-3p (Table S11). In our previous studies, CpFUL-like was regarded as a candidate gene for the regulation of wintersweet floral bud dormancy [41], and 35S::CpFUL-Like/Col-0 transgenic Arabidopsis lines showed an early flowering phenotype compared with WT (unpublished data).
The SVP/DAM gene plays a central regulatory role in bud dormancy [21][22][23][24][25][26][27]. In white pear, miRNA6390 was found to directly target PpDAM and regulate the release of floral bud dormancy [21]. In our previous study, two SVP homologous genes (CpSVP1 and CpSVP2) were identified from the transcriptome of wintersweet floral buds and considered as candidate genes for regulating floral bud dormancy [41]. However, miRNA6390 was not identified in all of the sRNA libraries in this study, and miRNAs targeting CpSVP1 or CpSVP2 were not found, which is similar to the results reported in Japanese pear, peach, and apple [37,39,59]. The reason for such differences may be different species, indicating that the regulatory roles of miRNAs in plant bud dormancy release may vary among species or even cultivars.
Plant hormones have been proved to be the most important internal factors regulating the dormancy cycle of deciduous trees [29]. Among them, the roles of ABA in dormancy regulation were most intensively studied, and numerous studies have shown that ABA plays a vital role in the establishment of bud dormancy [62][63][64][65][66]. In this study, of the 1707 target genes, some genes were involved in ABA biosynthesis and signal transduction, such as zeaxanthin epoxidase (ZEP), protein phosphatase 2c (PP2C), and abscisic acid-insensitive 5 (ABI5) ( Table S10). ZEP is an important regulatory enzyme in the first step of ABA biosynthesis, and the content of endogenous ABA is closely related to bud dormancy process [29]. PP2C, considered to be the inhibitor, is a vital component of the ABA signal transduction pathway [29]. In pear, PP2C is expressed at a low level during bud dormancy and upregulated after the release of bud dormancy [65,67]. ABI5 mediates ABA signal transduction. In our previous study, ABI5 (Unigene0029724) was regarded as a candidate gene for the regulation of floral bud dormancy in wintersweet. Its expression was significantly downregulated with the accumulation of low temperatures [41], and it was targeted by novel-m0128-3p (Table S11). Thus, it is speculated that miRNAs mediate ABA biosynthesis and signal transduction to regulate the dormancy release of floral buds in wintersweet.
Auxin, a major growth promoter, is involved in the regulation of dormancy release in many species. In tea plants, the content of endogenous IAA remains low during dormancy and increases continuously after dormancy release until spring germination [68]. In Prunus mume, the endogenous IAA level decreases as the floral buds enter dormancy and then increase continuously with the release of dormancy [69]. Furthermore, the transcriptome data of Prunus mume, Japanese apricot, and poplar show that the expression levels of auxin-associated genes are significantly different in buds between dormancy and dormancy release, and most of them are downregulated during dormancy [69][70][71]. However, for bud dormancy, the mechanism of auxin remains unclear. In this study, of the 1707 target genes, some genes involved in auxin signal transduction were found, including ARF, IAA27, AUX22, SAUR36, and TIR1 (Table S10). The auxin/indole-3-acetic acid (Aux/IAA) family, auxin response factor (ARF) family, and small auxin upregulated RNA (SAUR) are auxin response factors. Transport inhibitor response 1 (TIR1) is an auxin receptor and mediates transcriptional responses to auxin and the ubiquitination degradation of Aux/IAA [72][73][74]. The transcription factor ARF regulates the expression of auxin responsive genes and is regulated negatively by Aux/IAA [72,75]. These results indicate that miRNAs may regulate the dormancy release of floral buds by mediating auxin signal transduction in wintersweet.
Epigenetic regulation, such as DNA methylation and histone modification, is closely related to the induction, maintenance, and release of bud dormancy [18,[30][31][32][33][34][35]. Here, the 1707 target genes contained multiple DNA methyltransferases and demethylases, and they also contained many enzymes involved in histone modifications such as histone methylation, acetylation, and ubiquitination (Table S10), indicating that epigenetic modifications play important roles in the dormancy release of wintersweet floral buds. Among these genes, a histone acetyltransferase HAT-B and a histone deacetylase HDA9 were predicted to be targeted by miR2092-x and novel-m0128-3p, respectively (Table S11). Studies have shown that histone acetylation and deacetylation, which are catalyzed by histone acetyl-transferases (HATs) and histone deacetylases (HDAs), respectively, play regulatory roles in cold responses in plants [76].
In addition, in this study, miR157-y and miR6293-x were found to target RNA-dependent RNA polymerase 6 (RDR6) and Dicer-like 4 (DCL4), respectively (Table S11). RDR6 and DCL4 are mainly involved in the biogenesis of trans-acting siRNAs (ta-siRNAs) [77][78][79]. ta-siRNAs are a class of endogenous siRNAs with the length of 21 nt in plants, and they function similarly to miRNAs as post-transcriptional negative regulators [79,80]. It was reported that siRNAs could be involved in regulating the expression of MADS-box genes during bud dormancy transition in sweet cherry [35].
In this study, there were also a mass of target genes identified with unknown functions (Table S10). Furthermore, degradome sequencing could not obtain the information that miRNAs regulate target genes by translational repression. Therefore, more studies will be needed to clarify the regulation of miRNAs in wintersweet floral bud dormancy.

Plant Materials and Total RNA Extraction
All plant materials used in this study were floral buds (FBs) of 18-year-old potted C. praecox 'Concolor' plants that were grown on the campus of Southwest University (106 • 43 E, 29 • 83 N, Beibei District, Chongqing City, China) with different chilling accumulation. Chill units (CU) were calculated according to the "UTAH Model" [81]. Plant treatment and sample collection were as described in our previous study (Materials and methods, pages 2-3) [41]. Samples included FB.Nov (FBs in November before treatment), nF16 (non-flowering FBs collected after treatment at 16 • C for −300 CU), FB150, FB300, FB450, and IB570 (FBs initiate blooming) (FBs collected after treatment at 12 • C for 150, 300, 450, and 570 CU, respectively). Total RNA was extracted from each sample using Trizol reagent (Invitrogen, Waltham, MA, USA) by following the manufacturer's instructions.

sRNA Library Construction and Sequencing
After extraction of total RNA, RNA fragments of lengths 18-30 nt were enriched through polyacrylamide gel electrophoresis (PAGE) and were ligated with 3 and 5 adapters. Thereafter, reverse transcription and PCR amplification were performed, and 140-160 bp DNA fragments were enriched through agarose gel electrophoresis to generate the sRNA library [82]. The quality of the sRNA library was assessed using Agilent 2100 Bioanalyzer and qRT-PCR, and sequencing was performed using the Illumina HiSeq TM 2500 platform by Gene Denovo Biotechnology Co. (Guangzhou, China).

Degradome Library Construction and Sequencing
RNA fragments with 3 poly(A) were enriched from total RNA samples using magnetic beads with Oligo-d(T) and ligated with 5 adapters. Then, the first strand of cDNA was reverse-transcribed using biotinylated random primers. Subsequently, PCR amplification was performed to generate the library [11,[83][84][85]. Degradome sequencing was performed using the Illumina HiSeq 2000 platform by Gene Denovo Biotechnology Co. (Guangzhou, China) with the SE51 sequencing strategy to generate 50-nt reads.

Analysis of sRNA Sequencing Data
Raw reads from sequencing were filtered (removing contaminants, low-quality reads, and reads with length <18 nt reads) to obtain clean tags. All clean tags were mapped to the Genbank (Release 209.0, http://www.ncbi.nlm.nih.gov/genbank, accessed on 11 May 2019) [86], Rfam 11.0 (http://rfam.sanger.ac.uk/, accessed on 11 May 2019) [87], and RepeatMasker [88] databases to identify and remove non-coding structural RNAs (including rRNA, tRNA, snRNA, snoRNA, and scRNA) and possible repeat sequences. Subsequently, the clean tags were aligned with the C. praecox transcriptome database (NCBI accession number: PRJNA613935) established in our previous study [41] to remove the fragments from mRNA degradation. The remaining unannotated clean tags were used for miRNA identification.
The unannotated clean tags were aligned with miRBase 21.0 (http://www.mirbase.org, accessed on 11 May 2019) to identify known miRNAs in wintersweet samples with no more than two mismatches [89]. Thereafter, the remaining clean tags were mapped to the C. praecox transcriptome database to identify potentially novel miRNAs. Using default parameters, Mireap_v0.2 (http://sourceforge.net/projects/mireap/, accessed on 11 May 2019) was used to predict novel miRNAs and their hairpin structures [90,91]. The remaining tags were marked "unannotated".

Analysis of Differentially Expressed miRNAs
Based on the abundance of miRNAs in each sample, the miRNA expression level was calculated and normalized to transcripts per million (TPM) with the following formula: TPM = actual count of miRNA tags/total count of clean tags × 10 6 [11]. Following normalization, if the expression level of a given miRNA was zero, then the expression value was adjusted to 0.01 for further calculation [44,90]. Differentially expressed miRNAs (DEMs) were identified with edgeR (using default parameters) [92] with the criteria of |log 2 (FC)| ≥ 1 and p-value < 0.05. Fold change (FC) represents the ratio of the expression level between two samples. FC ≥ 2 and FC ≤ 0.5 represent upregulated and downregulated expression, respectively.

Analysis of Degradome Sequencing Data and Identification of miRNA Targets
After degradome sequencing was complete, adapters, contaminants, and low-quality reads were removed from raw reads to obtain clean tags. Then, non-coding structural RNAs were removed from clean tags as done for the sRNA-seq data. Degraded fragments with a single base ratio > 0.7 were defined as polyN, which were also removed from clean tags. The remaining clean tags were aligned to the C. praecox transcriptome database using the SOAP 2.20 (http://soap.genomics.org.cn/, accessed on 28 May 2019) [93] software to detect potentially degraded target genes. The tags mapped to the transcriptome were annotated as cDNA (cDNA_sense and cDNA_antisense) and used for the prediction of miRNA cleavage sites. CleaveLand 3.0 (http://axtell-lab-psu.weebly.com/cleaveland.html, accessed on 28 May 2019) [94] was used to count cleavage sites and build degradome density files. PAREsnip (http://srna-workbench.cmp.uea.ac.uk/tools/paresnip/, accessed on 28 May 2019) [95] was used to detect miRNA-target pairs with the criterion of p-value < 0.05. Finally, the target plot (T-plot) of miRNA-target pairs was constructed based on the abundance and position of degradome density files with CleaveLand 3.0. The process of degradome sequencing is shown in Figure S3. In addition, Gene Ontology (GO) enrichment analysis was performed for the identified target genes of miRNAs. The target genes were mapped to the GO database (http://www.geneontology.org/, accessed on 28 May 2019) using the GOseq R package, and the number of genes for each term was calculated. Then, the obtained statistical results were subjected to hypergeometric testing, and GO terms with p-value ≤ 0.05 were considered significantly enriched [96].
miRNA-target pairs were grouped into five categories according to the abundance of tags at the cleavage site with decreasing reliability [97]: category 0, the abundance of tags at this cleavage site is only the maximum on the transcript; category 1, the abundance of tags at this cleavage site and others is tied for the maximum on the transcript; category 2, the abundance of tags at this cleavage site is higher than the median but not the maximum on the transcript; category 3, the abundance of tags at this cleavage site is equal to or less than the median on the transcript; category 4, only one tag is found at this cleavage site.

Conclusions
This is the first study to analyze and identify the miRNAs and their target genes that regulate the dormancy release of floral buds in wintersweet using high-throughput sequencing technology. As a result, a total of 862 known and 402 novel miRNAs were identified with sRNA sequencing. The length of their sequences ranged from 18 nt to 26 nt, with 21-nt miRNAs being the most abundant (85, 50%), and the nucleotide at the 5 terminal of these miRNA was mainly biased toward U or A. Moreover, 23 DEMs (10 known and 13 novel) were screened by comparing the breaking and other dormancy samples, which indicated that miRNAs are involved in regulating the dormancy release of wintersweet floral buds. In addition, 1707 target genes of 21 DEMs were identified with degradome sequencing, and the potential functions of these targets were also annotated and analyzed preliminarily. The results provide valuable reference material for further elucidating the regulation of miRNAs in the dormancy release of wintersweet floral buds.